//========================================================================================
// Athena++ astrophysical MHD code
// Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
// Licensed under the 3-clause BSD License, see LICENSE file for details
//========================================================================================
//! \file gl_quadrature.cpp
//! \brief namespace containing GL quadrature function implementations, weights, and
//!  abscissas

// C headers

// C++ headers

// Athena++ headers
#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "gl_quadrature.hpp"

// Not checking that input integration order "N" is correct: 1 < N <=64

namespace GaussLegendre {
// NOTE: 3x function overload definitions follow the 64x abscissae and 64x weights arrays
// (around L1120)

// For n=2,4,16,20,32,40,48, refer to:
// Davis, P., & Rabinowitz, P. (1956). Abscissas and Weights for Gaussian Quadratures of
// High Order. Journal of Research of the National Bureau of Standards, 56(1).
static constexpr Real weight_n0[] = {0.0};
static constexpr Real weight_n1[] = {0.0};
static constexpr Real weight_n2[] = {
  1.0000000000000000000000,
};
static constexpr Real weight_n3[] = {
  0.8888888888888888888888,
  0.5555555555555555555555,
};
static constexpr Real weight_n4[] = {
  0.6521451548625461426269,
  0.3478548451374538573730,
};
static constexpr Real weight_n5[] = {
  0.5688888888888888888888,
  0.4786286704993664680412,
  0.2369268850561890875142,
};
static constexpr Real weight_n6[] = {
  0.3607615730481386075698,
  0.4679139345726910473898,
  0.1713244923791703450402,
};
static constexpr Real weight_n7[] = {
  0.4179591836734693877551,
  0.3818300505051189449503,
  0.2797053914892766679014,
  0.1294849661688696932706,
};
static constexpr Real weight_n8[] = {
  0.3626837833783619829651,
  0.3137066458778872873379,
  0.2223810344533744705443,
  0.1012285362903762591525,
};
static constexpr Real weight_n9[] = {
  0.3302393550012597631645,
  0.1806481606948574040584,
  0.0812743883615744119718,
  0.3123470770400028400686,
  0.2606106964029354623187,
};
static constexpr Real weight_n10[] = {
  0.2955242247147528701738,
  0.2692667193099963550912,
  0.2190863625159820439955,
  0.1494513491505805931457,
  0.0666713443086881375935,
};
static constexpr Real weight_n11[] = {
  0.2729250867779006307144,
  0.2628045445102466621806,
  0.2331937645919904799185,
  0.1862902109277342514260,
  0.1255803694649046246346,
  0.0556685671161736664827,
};
static constexpr Real weight_n12[] = {
  0.2491470458134027850005,
  0.2334925365383548087608,
  0.2031674267230659217490,
  0.1600783285433462263346,
  0.1069393259953184309602,
  0.0471753363865118271946,
};
static constexpr Real weight_n13[] = {
  0.2325515532308739101945,
  0.2262831802628972384120,
  0.2078160475368885023125,
  0.1781459807619457382800,
  0.1388735102197872384636,
  0.0921214998377284479144,
  0.0404840047653158795200,
};
static constexpr Real weight_n14[] = {
  0.2152638534631577901958,
  0.2051984637212956039659,
  0.1855383974779378137417,
  0.1572031671581935345696,
  0.1215185706879031846894,
  0.0801580871597602098056,
  0.0351194603317518630318,
};
static constexpr Real weight_n15[] = {
  0.2025782419255612728806,
  0.1984314853271115764561,
  0.1861610000155622110268,
  0.1662692058169939335532,
  0.1395706779261543144478,
  0.1071592204671719350118,
  0.0703660474881081247092,
  0.0307532419961172683546,
};
static constexpr Real weight_n16[] = {
  0.1894506104550684962853,
  0.1826034150449235888667,
  0.1691565193950025381893,
  0.1495959888165767320815,
  0.1246289712555338720524,
  0.0951585116824927848099,
  0.0622535239386478928628,
  0.0271524594117540948517,
};
static constexpr Real weight_n17[] = {
  0.1794464703562065254582,
  0.1765627053669926463252,
  0.1680041021564500445099,
  0.1540457610768102880814,
  0.1351363684685254732863,
  0.1118838471934039710947,
  0.0850361483171791808835,
  0.0554595293739872011294,
  0.0241483028685479319601,
};
static constexpr Real weight_n18[] = {
  0.1691423829631435918406,
  0.1642764837458327229860,
  0.1546846751262652449254,
  0.1406429146706506512047,
  0.1225552067114784601845,
  0.1009420441062871655628,
  0.0764257302548890565291,
  0.0497145488949697964533,
  0.0216160135264833103133,
};
static constexpr Real weight_n19[] = {
  0.1610544498487836959791,
  0.1589688433939543476499,
  0.1527660420658596667788,
  0.1426067021736066117757,
  0.1287539625393362276755,
  0.1115666455473339947160,
  0.0914900216224499994644,
  0.0690445427376412265807,
  0.0448142267656996003328,
  0.0194617882297264770363,
};
static constexpr Real weight_n20[] = {
  0.1527533871307258506980,
  0.1491729864726037467878,
  0.1420961093183820513292,
  0.1316886384491766268984,
  0.1181945319615184173123,
  0.1019301198172404350367,
  0.0832767415767047487247,
  0.0626720483341090635695,
  0.0406014298003869413310,
  0.0176140071391521183118,
};
static constexpr Real weight_n21[] = {
  0.1460811336496904271919,
  0.1445244039899700590638,
  0.1398873947910731547221,
  0.1322689386333374617810,
  0.1218314160537285341953,
  0.1087972991671483776634,
  0.0934444234560338615532,
  0.0761001136283793020170,
  0.0571344254268572082836,
  0.0369537897708524937999,
  0.0160172282577743333242,
};
static constexpr Real weight_n22[] = {
  0.1392518728556319933754,
  0.1365414983460151713525,
  0.1311735047870623707329,
  0.1232523768105124242855,
  0.1129322960805392183934,
  0.1004141444428809649320,
  0.0859416062170677274144,
  0.0697964684245204880949,
  0.0522933351526832859403,
  0.0337749015848141547933,
  0.0146279952982722006849,
};
static constexpr Real weight_n23[] = {
  0.1336545721861061753514,
  0.1324620394046966173716,
  0.1289057221880821499785,
  0.1230490843067295304675,
  0.1149966402224113649416,
  0.1048920914645414100740,
  0.0929157660600351474770,
  0.0792814117767189549228,
  0.0642324214085258521271,
  0.0480376717310846685716,
  0.0309880058569794443106,
  0.0134118594871417720813,
};
static constexpr Real weight_n24[] = {
  0.1279381953467521569740,
  0.1258374563468282961213,
  0.1216704729278033912044,
  0.1155056680537256013533,
  0.1074442701159656347825,
  0.0976186521041138882698,
  0.0861901615319532759171,
  0.0733464814110803057340,
  0.0592985849154367807463,
  0.0442774388174198061686,
  0.0285313886289336631813,
  0.0123412297999871995468,
};
static constexpr Real weight_n25[] = {
  0.1231760537267154512039,
  0.1222424429903100416889,
  0.1194557635357847722281,
  0.1148582591457116483393,
  0.1085196244742636531160,
  0.1005359490670506442022,
  0.0910282619829636498114,
  0.0801407003350010180132,
  0.0680383338123569172071,
  0.0549046959758351919259,
  0.0409391567013063126556,
  0.0263549866150321372619,
  0.0113937985010262879479,
};
static constexpr Real weight_n26[] = {
  0.1183214152792622765163, 0.1166604434852965820446,
  0.1133618165463196665494, 0.1084718405285765906565, 0.1020591610944254232384,
  0.0942138003559141484636, 0.0850458943134852392104, 0.0746841497656597458870,
  0.0632740463295748355394, 0.0509758252971478119983, 0.0379623832943627639503,
  0.0244178510926319087896, 0.0105513726173430071556,
};
static constexpr Real weight_n27[] = {
  0.1142208673789569890450, 0.1134763461089651486203, 0.1112524883568451926721,
  0.1075782857885331872121, 0.1025016378177457986712, 0.0960887273700285075656,
  0.0884231585437569501943, 0.0796048677730577712630, 0.0697488237662455929843,
  0.0589835368598335991103, 0.0474494125206150627040, 0.0352970537574197110225,
  0.0226862315961806231960, 0.0097989960512943602611,
};
static constexpr Real weight_n28[] = {
  0.1100470130164751962823, 0.1087111922582941352535, 0.1060557659228464179104,
  0.1021129675780607698142, 0.0969306579979299158504, 0.0905717443930328409421,
  0.0831134172289012183903, 0.0746462142345687790239, 0.0652729239669995957933,
  0.0551073456757167454314, 0.0442729347590042278395, 0.0329014277823043799776,
  0.0211321125927712597515, 0.0091242825930945177388,
};
static constexpr Real weight_n29[] = {
  0.1064793817183142442465, 0.1058761550973209414065, 0.1040733100777293739133,
  0.1010912737599149661218, 0.0969638340944086063019, 0.0917377571392587633479,
  0.0854722573661725275453, 0.0782383271357637838281, 0.0701179332550512785695,
  0.0612030906570791385421, 0.0515948269024979239125, 0.0414020625186828361048,
  0.0307404922020936226444, 0.0197320850561227059838, 0.0085169038787464096542,
};
static constexpr Real weight_n30[] = {
  0.1028526528935588403412, 0.1017623897484055045964, 0.0995934205867952670627,
  0.0963687371746442596394, 0.0921225222377861287176, 0.0868997872010829798023,
  0.0807558952294202153546, 0.0737559747377052062682, 0.0659742298821804951281,
  0.0574931562176190664817, 0.0484026728305940529029, 0.0387991925696270495968,
  0.0287847078833233693497, 0.0184664683110909591423, 0.0079681924961666056154,
};
static constexpr Real weight_n31[] = {
  0.0997205447934264514275, 0.0992250112266723078748, 0.0977433353863287250934,
  0.0952902429123195128072, 0.0918901138936414782153, 0.0875767406084778761261,
  0.0823929917615892639038, 0.0763903865987766164263, 0.0696285832354103661677,
  0.0621747865610284269103, 0.0541030824249168537116, 0.0454937075272011029023,
  0.0364322739123854640243, 0.0270090191849794218006, 0.0173186207903105824631,
  0.0074708315792487758586,
};
static constexpr Real weight_n32[] = {
  0.0965400885147278005667, 0.0956387200792748594190, 0.0938443990808045656391,
  0.0911738786957638847128, 0.0876520930044038111427, 0.0833119242269467552221,
  0.0781938957870703064717, 0.0723457941088485062253, 0.0658222227763618468376,
  0.0586840934785355471452, 0.0509980592623761761961, 0.0428358980222266806568,
  0.0342738629130214331026, 0.0253920653092620594557, 0.0162743947309056706051,
  0.0070186100094700966004,
};
static constexpr Real weight_n33[] = {
  0.0937684461602099965673, 0.0933564260655961161609, 0.0921239866433168462132,
  0.0900819586606385772397, 0.0872482876188443376072, 0.0836478760670387076139,
  0.0793123647948867383639, 0.0742798548439541493424, 0.0685945728186567128059,
  0.0623064825303174800316, 0.0554708466316635612849, 0.0481477428187116956701,
  0.0404015413316695915634, 0.0323003586323289532815, 0.0239155481017494803505,
  0.0153217015129346761279, 0.0066062278475873780586,
};
static constexpr Real weight_n34[] = {
  0.0909567403302598736153, 0.0902030443706407295739, 0.0887018978356938692870,
  0.0864657397470357497842, 0.0835130996998456551870, 0.0798684443397718447388,
  0.0755619746600319312708, 0.0706293758142557249990, 0.0651115215540764113785,
  0.0590541358275244931939, 0.0525074145726781061682, 0.0455256115233532724538,
  0.0381665937963875163217, 0.0304913806384461318094, 0.0225637219854949700840,
  0.0144501627485950354152, 0.0062291405559086847186,
};
static constexpr Real weight_n35[] = {
  0.0884867949071042906382, 0.0881405304302754629707, 0.0871044469971835342433,
  0.0853866533920991252259, 0.0830005937288565883799, 0.0799649422423242629326,
  0.0763034571554420535386, 0.0720447947725600646654, 0.0672222852690869039643,
  0.0618736719660801888870, 0.0560408162123701285783, 0.0497693704013535298051,
  0.0431084223261702187823, 0.0361101158634633805327, 0.0288292601088942540487,
  0.0213229799114835808834, 0.0136508283483614922664, 0.0058834334204430849757,
};
static constexpr Real weight_n36[] = {
  0.0859832756703947474900, 0.0853466857393386274918, 0.0840782189796619349334,
  0.0821872667043397095172, 0.0796878289120716019087, 0.0765984106458706745287,
  0.0729418850056530613538, 0.0687453238357364426136, 0.0640397973550154895563,
  0.0588601442453248173096, 0.0532447139777599190920, 0.0472350834902659784166,
  0.0408757509236448954741, 0.0342138107703072299212, 0.0272986214985687790944,
  0.0201815152977354715320, 0.0129159472840655744045, 0.0055657196642450453612,
};
static constexpr Real weight_n37[] = {
  0.0837683609931389047970, 0.0834745736258627872522, 0.0825952722364372508912,
  0.0811366245084650305098, 0.0791088618375293807672, 0.0765262075705292378858,
  0.0734067772484881727246, 0.0697724515557003448850, 0.0656487228727512494840,
  0.0610645165232259861309, 0.0560519879982749178085, 0.0506462976548246016038,
  0.0448853646624371666574, 0.0388096025019345444889, 0.0324616398475214810672,
  0.0258860369905589335227, 0.0191290444890839660435, 0.0122387801003075565263,
  0.0052730572794979393517,
};
static constexpr Real weight_n38[] = {
  0.0815250292803857866992, 0.0809824937705971006232, 0.0799010332435278215860,
  0.0782878446582109480753, 0.0761536635484463960659, 0.0735126925847434571452,
  0.0703825070668989547392, 0.0667839379791404119350, 0.0627409333921330540529,
  0.0582803991469972060223, 0.0534320199103323199737, 0.0482280618607586833743,
  0.0427031585046744342358, 0.0368940815940247381649, 0.0308395005451750546587,
  0.0245797397382323758952, 0.0181565777096132368988, 0.0116134447164686741776,
  0.0050028807496393456758,
};
static constexpr Real weight_n39[] = {
  0.0795276221394428524174, 0.0792762225683684710101, 0.0785236132873711767250,
  0.0772745525446820167285, 0.0755369373228360577047, 0.0733217534142686173811,
  0.0706430059706087607701, 0.0675176309662312653630, 0.0639653881386823889867,
  0.0600087360885961495749, 0.0556726903409162999073, 0.0509846652921294052140,
  0.0459743011089166318841, 0.0406732768479338439390, 0.0351151114981313307610,
  0.0293349559839033785921, 0.0233693848321781645947, 0.0172562290937249190408,
  0.0110347889391645942426, 0.0047529446916351013707,
};
static constexpr Real weight_n40[] = {
  0.0775059479784248112637, 0.0770398181642479655883, 0.0761103619006262423715,
  0.0747231690579682642001, 0.0728865823958040590605, 0.0706116473912867796954,
  0.0679120458152339038256, 0.0648040134566010380745, 0.0613062424929289391665,
  0.0574397690993915513666, 0.0532278469839368243549, 0.0486958076350722320614,
  0.0438709081856732719916, 0.0387821679744720176399, 0.0334601952825478473926,
  0.0279370069800234010984, 0.0222458491941669572615, 0.0164210583819078887128,
  0.0104982845311528136147, 0.0045212770985331912584,
};
static constexpr Real weight_n41[] = {
  0.0756955356472983723187, 0.0754787470927158240272, 0.0748296231762215518913,
  0.0737518820272234699392, 0.0722516968610230733963, 0.0703376606208174974816,
  0.0680207367608767667355, 0.0653141964535274104361, 0.0622335425809663164715,
  0.0587964209498719449911, 0.0550225192425787418801, 0.0509334542946174947811,
  0.0465526483690143420607, 0.0419051951959096894293, 0.0370177167035079884352,
  0.0319182117316992817870, 0.0266358992071104454675, 0.0212010633687795530756,
  0.0156449384078185885308, 0.0099999387739059453384, 0.0043061403581648876840,
};
static constexpr Real weight_n42[] = {
  0.0738642342321728799963, 0.0734608134534675282640, 0.0726561752438041048879,
  0.0714547142651709829218, 0.0698629924925941597661, 0.0678897033765219448553,
  0.0655456243649089789270, 0.0628435580450025764093, 0.0597982622275866543128,
  0.0564263693580183816464, 0.0527462956991740703439, 0.0487781407928032450274,
  0.0445435777719658778743, 0.0400657351806922617605, 0.0353690710975921108326,
  0.0304792406996034683629, 0.0254229595261130478867, 0.0202278695690526447570,
  0.0149224436973574941446, 0.0095362203017485024118, 0.0041059986046490846106,
};
static constexpr Real weight_n43[] = {
  0.0722157516937989879774, 0.0720275019714219743453, 0.0714637342525141412975,
  0.0705273877650850281262, 0.0692233441936566842822, 0.0675584022293651691924,
  0.0655412421263227974912, 0.0631823804493961123256, 0.0604941152499912945196,
  0.0574904619569105194276, 0.0541870803188817868633, 0.0506011927843901565238,
  0.0467514947543465800106, 0.0426580571979820837638, 0.0383422221941326575721,
  0.0338264920868602923449, 0.0291344132614984949159, 0.0242904566138388159020,
  0.0193199014236839003961, 0.0142487564315764861085, 0.0091039966374014033188,
  0.0039194902538441272829,
};
static constexpr Real weight_n44[] = {
  0.0705491577893540688113, 0.0701976854735582125871, 0.0694964918615725780370,
  0.0684490702693666609854, 0.0670606389062936523957, 0.0653381148791814349842,
  0.0632900797332038549501, 0.0609267367015619680385, 0.0582598598775954953342,
  0.0553027355637280525487, 0.0520700960917044618812, 0.0485780464483520375276,
  0.0448439840819700314462, 0.0408865123103462189084, 0.0367253478138088736429,
  0.0323812228120698208808, 0.0278757828212810100811, 0.0232314819020192106289,
  0.0184714817368147491720, 0.0136195867555799855202, 0.0087004813675248441225,
  0.0037454048031127775151,
};
static constexpr Real weight_n45[] = {
  0.0690418248292320201107, 0.0688773169776613228820, 0.0683845773786696745316,
  0.0675659541636075362709, 0.0664253484498425280829, 0.0649681957507234308538,
  0.0632014400738199377499, 0.0611335008310665225018, 0.0587742327188417385743,
  0.0561348787597864766439, 0.0532280167312689519459, 0.0500674992379520297991,
  0.0466683877183733652677, 0.0430468807091649711516, 0.0392202367293024475641,
  0.0352066922016090162476, 0.0310253749345154671625, 0.0266962139675776648056,
  0.0222398475505787323939, 0.0176775352579375906170, 0.0130311049915827843206,
  0.0083231892962182416457, 0.0035826631552835589311,
};
static constexpr Real weight_n46[] = {
  0.0675186858490364588202, 0.0672106136006781758623, 0.0665958747684548873757,
  0.0656772742677812073787, 0.0644590034671390695882, 0.0629466210643945081789,
  0.0611470277246504810153, 0.0590684345955463148075, 0.0567203258439912358168,
  0.0541134153858567544916, 0.0512595980071430213353, 0.0481718951017122005304,
  0.0448643952773181267670, 0.0413521901096787297042, 0.0376513053573860713276,
  0.0337786279991068965206, 0.0297518295522027557990, 0.0255892863971300106346,
  0.0213099987541365010544, 0.0169335140078362380462, 0.0124798837709886842067,
  0.0079698982297246224516, 0.0034303008681070482860,
};
static constexpr Real weight_n47[] = {
  0.0661351296236554796534, 0.0659905335888104745335, 0.0655573777665497402511,
  0.0648375562389457267026, 0.0638342166057170306312, 0.0625517462209216626405,
  0.0609957530087396453307, 0.0591730409423388759761, 0.0570915802932315402220,
  0.0547604727815302259571, 0.0521899117800571448722, 0.0493911377473611696045,
  0.0463763890865059112044, 0.0431588486484795382683, 0.0397525861225310037809,
  0.0361724965841749516134, 0.0324342355151847567676, 0.0285541507006433865047,
  0.0245492116596588185378, 0.0204369381476684276420, 0.0162353331464330596707,
  0.0119628484643123209639, 0.0076386162958488336141, 0.0032874538425280148832,
};
static constexpr Real weight_n48[] = {
  0.0647376968126839225030, 0.0644661644359500822065, 0.0639242385846481866239,
  0.0631141922862540256571, 0.0620394231598926639041, 0.0607044391658938800529,
  0.0591148396983956357464, 0.0572772921004032157051, 0.0551995036999841628682,
  0.0528901894851936670955, 0.0503590355538544749578, 0.0476166584924904748259,
  0.0446745608566942804194, 0.0415450829434647492140, 0.0382413510658307063172,
  0.0347772225647704388925, 0.0311672278327980889020, 0.0274265097083569482000,
  0.0235707608393243791405, 0.0196161604573555278144, 0.0155793157229438487281,
  0.0114772345792345394895, 0.0073275539012762621023, 0.0031533460523058386326,
};
static constexpr Real weight_n49[] = {
  0.0634632814047905977182, 0.0633355092964917485908, 0.0629527074651956994743,
  0.0623164173200572674010, 0.0614292009791929362968, 0.0602946309531520173031,
  0.0589172757600272660245, 0.0573026815301874754851, 0.0554573496748035886904,
  0.0533887107082589685279, 0.0511050943301445906746, 0.0486156958878282402776,
  0.0459305393555958535424, 0.0430604369812595979883, 0.0400169457663730213686,
  0.0368123209630006898194, 0.0334594667916221743424, 0.0299718846205838253506,
  0.0263636189270660169609, 0.0226492015874466764987, 0.0188435958530894584444,
  0.0149621449356246510295, 0.0110205510315935804975, 0.0070350995900864514734,
  0.0030272789889229050774,
};
static constexpr Real weight_n50[] = {
  0.0621766166553472623210, 0.0619360674206832433840, 0.0614558995903166637564,
  0.0607379708417702160317, 0.0597850587042654575095, 0.0586008498132224458351,
  0.0571899256477283837230, 0.0555577448062125176235, 0.0537106218889962465234,
  0.0516557030695811384899, 0.0494009384494663149212, 0.0469550513039484329656,
  0.0443275043388032754920, 0.0415284630901476974224, 0.0385687566125876752447,
  0.0354598356151461541607, 0.0322137282235780166481, 0.0288429935805351980299,
  0.0253606735700123904401, 0.0217802431701247929815, 0.0181155607134893903512,
  0.0143808227614855744193, 0.0105905483836509692635, 0.0067597991957454015027,
  0.0029086225531551409584,
};
static constexpr Real weight_n51[] = {
  0.0609989248412058801597, 0.0608854648448563438811, 0.0605455069347377951381,
  0.0599803157775032520900, 0.0591919939229615437835, 0.0581834739825921405984,
  0.0569585077202586621000, 0.0555216520957386930167, 0.0538782523130455614340,
  0.0520344219366970875641, 0.0499970201500574097795, 0.0477736262406231019999,
  0.0453725114076500687481, 0.0428026079978800866536, 0.0400734762854964531868,
  0.0371952689232602928429, 0.0341786932041883362362, 0.0310349712901600084544,
  0.0277757985941624771959, 0.0244133005737814342731, 0.0209599884017032105797,
  0.0174287147234010522595, 0.0138326340064778222966, 0.0101851912978217299392,
  0.0065003377832526002921, 0.0027968071710898955755,
};
static constexpr Real weight_n52[] = {
  0.0598103657452918602477, 0.0595962601712481582583, 0.0591688154660429703693,
  0.0585295617718138685502, 0.0576807874525268276539, 0.0566255309023685971908,
  0.0553675696693026525490, 0.0539114069327572647508, 0.0522622553839069930343,
  0.0504260185663423772182, 0.0484092697440748968539, 0.0462192283727847935076,
  0.0438637342590004079951, 0.0413512195005602716790, 0.0386906783104239789851,
  0.0358916348350972329419, 0.0329641090897187979150, 0.0299185811471439466412,
  0.0267659537465040134494, 0.0235175135539844615903, 0.0201848915079807922029,
  0.0167800233963007356779, 0.0133151149823409606566, 0.0098026345794627520619,
  0.0062555239629732768997, 0.0026913169500471111189,
};
static constexpr Real weight_n53[] = {
  0.0587187941511643645254, 0.0586175862327202633180, 0.0583143113622560075562,
  0.0578100149917131963196, 0.0571064355362671917733, 0.0562059983817397098086,
  0.0551118075239335990023, 0.0538276348687310290420, 0.0523579072298727181997,
  0.0507076910692927152964, 0.0488826750326991404204, 0.0468891503407503140218,
  0.0447339891036728102127, 0.0424246206345200135922, 0.0399690058435403821270,
  0.0373756098034829156741, 0.0346533725835342379583, 0.0318116784590193230632,
  0.0288603236178237362627, 0.0258094825107575177139, 0.0226696730570702083987,
  0.0194517211076368953880, 0.0161667252566874639280, 0.0128260261442403791791,
  0.0094412022849403443866, 0.0060242762269486732812, 0.0025916837205670318115,
};
static constexpr Real weight_n54[] = {
  0.0576175367071470246723, 0.0574261370541121148592, 0.0570439735587945985678,
  0.0564723157306259650310, 0.0557130625605899876833, 0.0547687362130579863062,
  0.0536424736475536112721, 0.0523380161982987446655, 0.0508596971461881443197,
  0.0492124273245288860687, 0.0474016788064449910585, 0.0454334667282767139748,
  0.0433143293095970154419, 0.0410513061366449742217, 0.0386519147821025168368,
  0.0361241258403835525828, 0.0334763364643726457160, 0.0307173424978706760540,
  0.0278563093105958702870, 0.0249027414672087730500, 0.0218664514228530859455,
  0.0187575276214693779120, 0.0155863030359241317029, 0.0123633281288476441664,
  0.0090993694555093969480, 0.0058056110152399848788, 0.0024974818357615857759,
};
static constexpr Real weight_n55[] = {
  0.0566029764445604254401, 0.0565123182497720014006, 0.0562406340710843680282,
  0.0557887941952840871029, 0.0551582460025086875966, 0.0543510093299111020703,
  0.0533696700016054727235, 0.0522173715456320845643, 0.0508978051244939792247,
  0.0494151977115517394807, 0.0477742985512006955500, 0.0459803639462838381039,
  0.0440391404216065898951, 0.0419568463177187623952, 0.0397401518743371796094,
  0.0373961578679655452829, 0.0349323728735898874072, 0.0323566892261858316847,
  0.0296773577651610412212, 0.0269029614563962706671, 0.0240423880097256220077,
  0.0211048016680164541202, 0.0180996145207290624079, 0.0150364583335117882131,
  0.0119251607198486121707, 0.0087757461070585281773, 0.0055986322665607673540,
  0.0024083236199797888191,
};
static constexpr Real weight_n56[] = {
  0.0555797463065143958462, 0.0554079525032451232177, 0.0550648959017624257963,
  0.0545516368708894210617, 0.0538697618657144857089, 0.0530213785240107639679,
  0.0520091091517413998430, 0.0508360826177984805601, 0.0495059246830475789199,
  0.0480227467936002581207, 0.0463911333730018967621, 0.0446161276526922832134,
  0.0427032160846670865110, 0.0406583113847445178801, 0.0384877342592476624868,
  0.0361981938723151860358, 0.0337967671156117612954, 0.0312908767473104478678,
  0.0286882684738227417298, 0.0259969870583919521918, 0.0232253515625653169372,
  0.0203819298824025726348, 0.0174755129114009465049, 0.0145150892780214718077,
  0.0115098243403833821737, 0.0084690631633078876616, 0.0054025222460153377613,
  0.0023238553757732155011,
};
static constexpr Real weight_n57[] = {
  0.0546343287565840240628, 0.0545528036047618864801, 0.0543084714524986431387,
  0.0539020614832985746428, 0.0533347865848191584265, 0.0526083397291774324402,
  0.0517248889205178247206, 0.0506870707249274086566, 0.0494979824020196789938,
  0.0481611726616877512688, 0.0466806310736415037838, 0.0450607761613811577972,
  0.0433064422162151965964, 0.0414228648708011103631, 0.0394156654754801140899,
  0.0372908343244173173547, 0.0350547127823126175057, 0.0327139743663715685424,
  0.0302756048426939994584, 0.0277468814021801923212, 0.0251353509909181226472,
  0.0224488078907764380796, 0.0196952706994885203824, 0.0168829590234415490350,
  0.0140202707907535561702, 0.0111157637323359901456, 0.0081781600678212326262,
  0.0052165334747187793905, 0.0022437538722506629097,
};
static constexpr Real weight_n58[] = {
  0.0536811198633348488639, 0.0535263433040582521006, 0.0532172364465790141034,
  0.0527546905263708334296, 0.0521400391836698189712, 0.0513750546182857254745,
  0.0504619424799531252976, 0.0494033355089623928665, 0.0482022859454177484065,
  0.0468622567290263469184, 0.0453871115148198025039, 0.0437811035336402510390,
  0.0420488633295821259945, 0.0401953854098677968880, 0.0382260138458584332294,
  0.0361464268670872705407, 0.0339626204934160107977, 0.0316808912538093273202,
  0.0293078180441604907183, 0.0268502431819818684759, 0.0243152527249639525402,
  0.0217101561401462357669, 0.0190424654618934086557, 0.0163198742349709650521,
  0.0135502371129888121451, 0.0107415535328787741168, 0.0079019738499986747540,
  0.0050399816126502430850, 0.0021677232496274499430,
};
static constexpr Real weight_n59[] = {
  0.0527980126219904214155, 0.0527244338591279319613, 0.0525039026478287390509,
  0.0521370336483753913839, 0.0516248493908914821464, 0.0509687774253939168502,
  0.0501706463429969028107, 0.0492326806793619857796, 0.0481574947146064403881,
  0.0469480851869620191931, 0.0456078229405097698318, 0.0441404435302973806907,
  0.0425500368110676386673, 0.0408410355386867076602, 0.0390182030161600095030,
  0.0370866198188709226918, 0.0350516696364001087837, 0.0329190242710452777575,
  0.0306946278361116832397, 0.0283846802005347979051, 0.0259956197312985001866,
  0.0235341053937133634252, 0.0210069982884371873504, 0.0184213427536100293606,
  0.0157843473130814661473, 0.0131033663063451910183, 0.0103858855009958621937,
  0.0076395294534875751427, 0.0048722391682652847685, 0.0020954922845412234027,
};
static constexpr Real weight_n60[] = {
  0.0519078776312206397328, 0.0517679431749101875438, 0.0514884515009809339950,
  0.0510701560698556274045, 0.0505141845325093745982, 0.0498220356905501810111,
  0.0489955754557568353894, 0.0480370318199711809636, 0.0469489888489122048470,
  0.0457343797161144866471, 0.0443964787957871133277, 0.0429388928359356419542,
  0.0413655512355847556131, 0.0396806954523807994701, 0.0378888675692434440309,
  0.0359948980510845030665, 0.0340038927249464228349, 0.0319212190192963289494,
  0.0297524915007889452408, 0.0275035567499247916352, 0.0251804776215212483795,
  0.0227895169439978198637, 0.0203371207294572867750, 0.0178299010142077202603,
  0.0152746185967847993067, 0.0126781664768159601314, 0.0100475571822879843578,
  0.0073899311633454555315, 0.0047127299269535686408, 0.0020268119688737584964,
};
static constexpr Real weight_n61[] = {
  0.0510811194407862179779, 0.0510144870386972635437, 0.0508147636688183432077,
  0.0504824703867974046481, 0.0500184741081782534250, 0.0494239853467355899399,
  0.0487005550564115260875, 0.0478500705850956071618, 0.0468747507508090659764,
  0.0457771400531459593713, 0.0445601020350834882715, 0.0432268118124960979010,
  0.0417807477908884920666, 0.0402256825909982473676, 0.0385656732070081727461,
  0.0368050504231548173843, 0.0349484075165333510908, 0.0330005882759074106327,
  0.0309666743683973948246, 0.0288519720881834015043, 0.0266619985241508896628,
  0.0244024671875442029153, 0.0220792731483190440024, 0.0196984777461011813305,
  0.0172662929876137435944, 0.0147890658849379145461, 0.0122732635078121046292,
  0.0097254618303561337361, 0.0071523549917490895858, 0.0045609240060124171845,
  0.0019614533616702826717,
};
static constexpr Real weight_n62[] = {
  0.0502480003752562816884, 0.0501210695690432880748, 0.0498675285949523942447,
  0.0494880179196992925278, 0.0489834962205178371048, 0.0483552379634776728348,
  0.0476048301841012322704, 0.0467341684784155248022, 0.0457454522145701807772,
  0.0446411789771244142936, 0.0434241382580474195800, 0.0420974044103850966430,
  0.0406643288824174409682, 0.0391285317519630841233, 0.0374938925822800299856,
  0.0357645406227681412855, 0.0339448443794105450911, 0.0320394005816246781063,
  0.0300530225739898700770, 0.0279907281633146375412, 0.0258577269540246980270,
  0.0236594072086827925745, 0.0214013222776699688411, 0.0190891766585731987325,
  0.0167288117901773162885, 0.0143261918238065177674, 0.0118873901170105019448,
  0.0094185794284203876379, 0.0069260419018309608717, 0.0044163334569309048132,
  0.0018992056795136904803,
};
static constexpr Real weight_n63[] = {
  0.0494723666239310208886, 0.0494118330399181789670, 0.0492303804237475607850,
  0.0489284528205119899447, 0.0485067890978838478640, 0.0479664211379951314110,
  0.0473086713122689190806, 0.0465351492453836965103, 0.0456477478762926086858,
  0.0446486388259413953703, 0.0435402670830275907989, 0.0423253450208158229825,
  0.0410068457596663986351, 0.0395879958915440939848, 0.0380722675843495567636,
  0.0364633700854572896304, 0.0347652406453558776971, 0.0329820348837793417656,
  0.0311181166222198175082, 0.0291780472082805269455, 0.0271665743590979332251,
  0.0250886205533449866186, 0.0229492710048899331489, 0.0207537612580390907753,
  0.0185074644601612704092, 0.0162158784103383388822, 0.0138846126161156108248,
  0.0115193760768800417507, 0.0091259686763266563540, 0.0067102917659601362519,
  0.0042785083468637618660, 0.0018398745955770841170,
};
static constexpr Real weight_n64[] = {
  0.0486909570091397203833, 0.0485754674415034269347, 0.0483447622348029571697,
  0.0479993885964583077281, 0.0475401657148303086622, 0.0469681828162100173253,
  0.0462847965813144172959, 0.0454916279274181444797, 0.0445905581637565630601,
  0.0435837245293234533768, 0.0424735151236535890073, 0.0412625632426235286101,
  0.0399537411327203413866, 0.0385501531786156291289, 0.0370551285402400460404,
  0.0354722132568823838106, 0.0338051618371416093915, 0.0320579283548515535854,
  0.0302346570724024788679, 0.0283396726142594832275, 0.0263774697150546586716,
  0.0243527025687108733381, 0.0222701738083832541592, 0.0201348231535302093723,
  0.0179517157756973430850, 0.0157260304760247193219, 0.0134630478967186425980,
  0.0111681394601311288185, 0.0088467598263639477230, 0.0065044579689783628561,
  0.0041470332605624676352, 0.0017832807216964329472,
};


static constexpr Real abscissa_n0[] = {-100.};
static constexpr Real abscissa_n1[] = {-100.};
static constexpr Real abscissa_n2[] = {
  0.5773502691896257645091
};
static constexpr Real abscissa_n3[] = {
  0.0, 0.7745966692414833770358
};
static constexpr Real abscissa_n4[] = {
  0.3399810435848562648026, 0.8611363115940525752239
};
static constexpr Real abscissa_n5[] = {
  0.0, 0.5384693101056830910363, 0.9061798459386639927976
};
static constexpr Real abscissa_n6[] = {
  0.6612093864662645136613, 0.2386191860831969086305, 0.9324695142031520278123
};
static constexpr Real abscissa_n7[] = {
  0.0, 0.4058451513773971669066, 0.7415311855993944398638, 0.9491079123427585245261
};
static constexpr Real abscissa_n8[] = {
  0.1834346424956498049394, 0.5255324099163289858177, 0.7966664774136267395915,
  0.9602898564975362316835
};
static constexpr Real abscissa_n9[] = {
  0.0, 0.8360311073266357942994, 0.9681602395076260898355, 0.3242534234038089290385,
  0.6133714327005903973087
};
static constexpr Real abscissa_n10[] = {
  0.1488743389816312108848, 0.4333953941292471907992, 0.6794095682990244062343,
  0.8650633666889845107320, 0.9739065285171717200779
};
static constexpr Real abscissa_n11[] = {
  0.0, 0.2695431559523449723315, 0.5190961292068118159257, 0.7301520055740493240934,
  0.8870625997680952990751, 0.9782286581460569928039
};
static constexpr Real abscissa_n12[] = {
  0.1252334085114689154724, 0.3678314989981801937526, 0.5873179542866174472967,
  0.7699026741943046870368, 0.9041172563704748566784, 0.9815606342467192506905
};
static constexpr Real abscissa_n13[] = {
  0.0, 0.2304583159551347940655, 0.4484927510364468528779, 0.6423493394403402206439,
  0.8015780907333099127942, 0.9175983992229779652065, 0.9841830547185881494728
};
static constexpr Real abscissa_n14[] = {
  0.1080549487073436620662, 0.3191123689278897604356, 0.5152486363581540919652,
  0.6872929048116854701480, 0.8272013150697649931897, 0.9284348836635735173363,
  0.9862838086968123388415
};
static constexpr Real abscissa_n15[] = {
  0.0, 0.2011940939974345223006, 0.3941513470775633698972, 0.5709721726085388475372,
  0.7244177313601700474161, 0.8482065834104272162006, 0.9372733924007059043077,
  0.9879925180204854284895
};
static constexpr Real abscissa_n16[] = {
  0.0950125098376374401853, 0.2816035507792589132304, 0.4580167776572273863424,
  0.6178762444026437484466, 0.7554044083550030338951, 0.8656312023878317438804,
  0.9445750230732325760779, 0.9894009349916499325961
};
static constexpr Real abscissa_n17[] = {
  0.0, 0.1784841814958478558506, 0.3512317634538763152971, 0.5126905370864769678862,
  0.6576711592166907658503, 0.7815140038968014069252, 0.8802391537269859021229,
  0.9506755217687677612227, 0.9905754753144173356754
};
static constexpr Real abscissa_n18[] = {
  0.0847750130417353012422, 0.2518862256915055095889, 0.4117511614628426460359,
  0.5597708310739475346078, 0.6916870430603532078748, 0.8037049589725231156824,
  0.8926024664975557392060, 0.9558239495713977551811, 0.9915651684209309467300
};
static constexpr Real abscissa_n19[] = {
  0.0, 0.1603586456402253758680, 0.3165640999636298319901, 0.4645707413759609457172,
  0.6005453046616810234696, 0.7209661773352293786170, 0.8227146565371428249789,
  0.9031559036148179016426, 0.9602081521348300308527, 0.9924068438435844031890
};
static constexpr Real abscissa_n20[] = {
  0.0765265211334973337546, 0.2277858511416450780804, 0.3737060887154195606725,
  0.5108670019508270980043, 0.6360536807265150254528, 0.7463319064601507926143,
  0.8391169718222188233945, 0.9122344282513259058677, 0.9639719272779137912676,
  0.9931285991850949247861
};
static constexpr Real abscissa_n21[] = {
  0.0, 0.1455618541608950909370, 0.2880213168024010966007, 0.4243421202074387835736,
  0.5516188358872198070590, 0.6671388041974123193059, 0.7684399634756779086158,
  0.8533633645833172836472, 0.9200993341504008287901, 0.9672268385663062943166,
  0.9937521706203895002602
};
static constexpr Real abscissa_n22[] = {
  0.0697392733197222212138, 0.2078604266882212854788, 0.3419358208920842251581,
  0.4693558379867570264063, 0.5876404035069115929588, 0.6944872631866827800506,
  0.7878168059792081620042, 0.8658125777203001365364, 0.9269567721871740005206,
  0.9700604978354287271239, 0.9942945854823992920730
};
static constexpr Real abscissa_n23[] = {
  0.0, 0.1332568242984661109317, 0.2641356809703449305338, 0.3903010380302908314214,
  0.5095014778460075496897, 0.6196098757636461563850, 0.7186613631319501944616,
  0.8048884016188398921511, 0.8767523582704416673781, 0.9329710868260161023491,
  0.9725424712181152319560, 0.9947693349975521235239
};
static constexpr Real abscissa_n24[] = {
  0.0640568928626056260850, 0.1911188674736163091586, 0.3150426796961633743867,
  0.4337935076260451384870, 0.5454214713888395356583, 0.6480936519369755692524,
  0.7401241915785543642438, 0.8200019859739029219539, 0.8864155270044010342131,
  0.9382745520027327585236, 0.9747285559713094981983, 0.9951872199970213601799
};
static constexpr Real abscissa_n25[] = {
  0.0, 0.1228646926107103963873, 0.2438668837209884320451, 0.3611723058093878377358,
  0.4730027314457149605221, 0.5776629302412229677236, 0.6735663684734683644851,
  0.7592592630373576305772, 0.8334426287608340014210, 0.8949919978782753688510,
  0.9429745712289743394140, 0.9766639214595175114983, 0.9955569697904980979087
};
static constexpr Real abscissa_n26[] = {
  0.0592300934293132070937, 0.1768588203568901839690, 0.2920048394859568951428,
  0.4030517551234863064810, 0.5084407148245057176957, 0.6066922930176180632319,
  0.6964272604199572648638, 0.7763859488206788561929, 0.8454459427884980187975,
  0.9026378619843070742176, 0.9471590666617142501359, 0.9783854459564709911005,
  0.9958857011456169290032
};
static constexpr Real abscissa_n27[] = {
  0.0, 0.1139725856095299669328, 0.2264593654395368588572, 0.3359939036385088997303,
  0.4411482517500268805859, 0.5405515645794568949003, 0.6329079719464951409277,
  0.7170134737394236992948, 0.7917716390705082271443, 0.8562079080182944903027,
  0.9094823206774911043006, 0.9509005578147050068519, 0.9799234759615012228558,
  0.9961792628889885669388
};
static constexpr Real abscissa_n28[] = {
  0.0550792898840342704265, 0.1645692821333807712814, 0.2720616276351780776768,
  0.3762515160890787102213, 0.4758742249551182610344, 0.5697204718114017193080,
  0.6566510940388649612198, 0.7356108780136317720281, 0.8056413709171791714478,
  0.8658925225743950489422, 0.9156330263921320738696, 0.9542592806289381972541,
  0.9813031653708727536945, 0.9964424975739544499504
};
static constexpr Real abscissa_n29[] = {
  0.0, 0.1062782301326792301709, 0.2113522861660010745063, 0.3140316378676399349481,
  0.4131528881740086638907, 0.5075929551242276421026, 0.5962817971382278203795,
  0.6782145376026865151561, 0.7524628517344771339126, 0.8181854876152524449895,
  0.8746378049201027904177, 0.9211802329530587850937, 0.9572855957780877257982,
  0.9825455052614131748709, 0.9966794422605965861631
};
static constexpr Real abscissa_n30[] = {
  0.0514718425553176958330, 0.1538699136085835469637, 0.2546369261678898464398,
  0.3527047255308781134710, 0.4470337695380891767806, 0.5366241481420198992641,
  0.6205261829892428611404, 0.6978504947933157969322, 0.7677774321048261949179,
  0.8295657623827683974428, 0.8825605357920526815431, 0.9262000474292743258793,
  0.9600218649683075122168, 0.9836681232797472099700, 0.9968934840746495402716
};
static constexpr Real abscissa_n31[] = {
  0.0, 0.0995553121523415203251, 0.1981211993355706287724, 0.2947180699817016166179,
  0.3883859016082329430613, 0.4781937820449024804405, 0.5632491614071492627209,
  0.6427067229242603461844, 0.7157767845868532839059, 0.7817331484166249404063,
  0.8399203201462673400869, 0.8897600299482710433741, 0.9307569978966481649569,
  0.9625039250929496617890, 0.9846859096651524840024, 0.9970874818194770740556
};
static constexpr Real abscissa_n32[] = {
  0.0483076656877383162348, 0.1444719615827964934851, 0.2392873622521370745446,
  0.3318686022821276497799, 0.4213512761306353453641, 0.5068999089322293900237,
  0.5877157572407623290407, 0.6630442669302152009751, 0.7321821187402896803874,
  0.7944837959679424069630, 0.8493676137325699701336, 0.8963211557660521239653,
  0.9349060759377396891709, 0.9647622555875064307738, 0.9856115115452683354001,
  0.9972638618494815635449
};
static constexpr Real abscissa_n33[] = {
  0.0, 0.0936310658547333856707, 0.1864392988279915723357, 0.2776090971524970294032,
  0.3663392577480733410702, 0.4518500172724506957259, 0.5333899047863476435488,
  0.6102423458363790273072, 0.6817319599697427862682, 0.7472304964495621578590,
  0.8061623562741665897962, 0.8580096526765040646430, 0.9023167677434335830405,
  0.9386943726111683503558, 0.9668229096899927689283, 0.9864557262306424881103,
  0.9974246942464552172661
};
static constexpr Real abscissa_n34[] = {
  0.0455098219531025427490, 0.1361523572591829758944, 0.2256666916164494838686,
  0.3133110813394632474583, 0.3983592777586459406314, 0.4801065451903270341941,
  0.5578755006697466427364, 0.6310217270805285453177, 0.6989391132162629079330,
  0.7610648766298730141874, 0.8168842279009336645915, 0.8659346383345644692635,
  0.9078096777183244688008, 0.9421623974051070916316, 0.9687082625333442817646,
  0.9872278164063094850497, 0.9975717537908419192433
};
static constexpr Real abscissa_n35[] = {
  0.0, 0.0883713432756592636009, 0.1760510611659895699743, 0.2623529412092960579708,
  0.3466015544308139458769, 0.4281375415178142541876, 0.5063227732414886150242,
  0.5805453447497645099345, 0.6502243646658903886757, 0.7148145015566287832644,
  0.7738102522869125552674, 0.8267498990922254068340, 0.8732191250252223315232,
  0.9128542613593176144649, 0.9453451482078273295387, 0.9704376160392298332150,
  0.9879357644438514980351, 0.9977065690996002972601
};
static constexpr Real abscissa_n36[] = {
  0.0430181984737086072269, 0.1287361038093847886519, 0.2135008923168655789432,
  0.2966849953440282705032, 0.3776725471196892163227, 0.4558639444334202672072,
  0.5306802859262451616407, 0.6015676581359805350799, 0.6680012365855210620971,
  0.7294891715935565820902, 0.7855762301322065128277, 0.8358471669924753064188,
  0.8799298008903971319824, 0.9174977745156590660758, 0.9482729843995075452024,
  0.9720276910496979493356, 0.9885864789022122380733, 0.9978304624840858361988
};
static constexpr Real abscissa_n37[] = {
  0.0, 0.0836704089547699019430, 0.1667539302398519769686, 0.2486677927913657588055,
  0.3288374298837069994984, 0.4067005093183261101009, 0.4817108778032055541471,
  0.5533423918615817812350, 0.6210926084089244831481, 0.6844863091309593574455,
  0.7430788339819652625470, 0.7964592005099022933930, 0.8442529873405559679874,
  0.8861249621554860789461, 0.9217814374124637426681, 0.9509723432620948213285,
  0.9734930300564857443288, 0.9891859632143191866840, 0.9979445824779136489408
};
static constexpr Real abscissa_n38[] = {
  0.0407851479045782399133, 0.1220840253378674198696, 0.2025704538921167032039,
  0.2817088097901652613601, 0.3589724404794350132567, 0.4338471694323764843732,
  0.5058347179279311032405, 0.5744560210478070811329, 0.6392544158296817071803,
  0.6997986803791843559128, 0.7556859037539706807377, 0.8065441676053168155515,
  0.8520350219323621888596, 0.8918557390046322167949, 0.9257413320485843968251,
  0.9534663309335295956705, 0.9748463285901535076408, 0.9897394542663855719444,
  0.9980499305356876198128
};
static constexpr Real abscissa_n39[] = {
  0.0, 0.0794438046087554775819, 0.1583853399978377999227, 0.2363255124618357673360,
  0.3127715592481859225359, 0.3872401639715614558538, 0.4592605123091360486632,
  0.5283772686604374738963, 0.5941534549572779886928, 0.6561732134320109107344,
  0.7140444358945346791338, 0.7674012429310634998322, 0.8159062974301431043532,
  0.8592529379999061539137, 0.8971671192929928878482, 0.9294091484867382296978,
  0.9557752123246522771108, 0.9760987093334710538448, 0.9902515368546859836397,
  0.9981473830664329060054
};
static constexpr Real abscissa_n40[] = {
  0.0387724175060508219331, 0.1160840706752552084834, 0.1926975807013710997155,
  0.2681521850072536811411, 0.3419940908257584730074, 0.4137792043716050015248,
  0.4830758016861787129085, 0.5494671250951282020759, 0.6125538896679802379526,
  0.6719566846141795483793, 0.7273182551899271032809, 0.7783056514265193876949,
  0.8246122308333116631963, 0.8659595032122595038207, 0.9020988069688742967282,
  0.9328128082786765333608, 0.9579168192137916558045, 0.9772599499837742626633,
  0.9907262386994570064530, 0.9982377097105592003496
};
static constexpr Real abscissa_n41[] = {
  0.0, 0.0756232589891629969237, 0.1508133548639921635743, 0.2251396056334227756057,
  0.2981762773418248659229, 0.3695050226404814414283, 0.4387172770514070885171,
  0.5054165991994060327083, 0.5692209416102158696547, 0.6297648390721963204886,
  0.6867015020349512895846, 0.7397048030699261810601, 0.7884711450474093727362,
  0.8327212004013613312442, 0.8722015116924414088336, 0.9066859447581011729583,
  0.9359769874978538256823, 0.9599068917303462260994, 0.9783386735610833844691,
  0.9911671096990163082501, 0.9983215885747714415191
};
static constexpr Real abscissa_n42[] = {
  0.0369489431653517758130, 0.1106450272085198683491, 0.1837368065648545508527,
  0.2558250793428790839664, 0.3265161244654115121971, 0.3954238520429750576770,
  0.4621719120704219297590, 0.5263957499311922875928, 0.5877445974851093228407,
  0.6458833888692478339574, 0.7004945905561712137415, 0.7512799356894804895684,
  0.7979620532554874132327, 0.8402859832618169009254, 0.8780205698121727427119,
  0.9109597249041274525838, 0.9389235573549881785331, 0.9617593653382044887469,
  0.9793425080637481937089, 0.9915772883408609197923, 0.9983996189900624150228
};
static constexpr Real abscissa_n43[] = {
  0.0, 0.0721529908745862354223, 0.1439298095107133107702, 0.2149562448605182090148,
  0.2848619980329136271057, 0.3532826128643038066454, 0.4198613760292692524872,
  0.4842511767857347240696, 0.5461163166600847191404, 0.6051342596396009357250,
  0.6609973137514981331651, 0.7134142352689570548519, 0.7621117471949551214603,
  0.8068359641369386352790, 0.8473537162093150489990, 0.8834537652186168633376,
  0.9149479072061387294561, 0.9416719568476378618184, 0.9634866130140799934101,
  0.9802782209802553315055, 0.9919595575932441464207, 0.9984723322425077135176
};
static constexpr Real abscissa_n44[] = {
  0.0352892369641353590581, 0.1056919017086532471173, 0.1755680147755167857465,
  0.2445694569282012515073, 0.3123524665027858122365, 0.3785793520147071325117,
  0.4429201745254114838348, 0.5050543913882023179827, 0.5646724531854707684246,
  0.6214773459035758478024, 0.6751860706661223653336, 0.7255310536607170026069,
  0.7722614792487558990177, 0.8151445396451350104874, 0.8539665950047103787283,
  0.8885342382860432023383, 0.9186752599841757743230, 0.9442395091181940992032,
  0.9650996504224931393943, 0.9811518330779139666627, 0.9923163921385158084833,
  0.9985402006367742249360
};
static constexpr Real abscissa_n45[] = {
  0.0, 0.0689869801631441724904, 0.1376452059832530287565, 0.2056474897832637457197,
  0.2726697697523775606087, 0.3383926542506021616434, 0.4025029438585419140779,
  0.4646951239196350985796, 0.5246728204629160670911, 0.5821502125693531866809,
  0.6368533944532233592712, 0.6885216807712005252320, 0.7369088489454903526237,
  0.7817843125939062913123, 0.8229342205020863370357, 0.8601624759606642253390,
  0.8932916717532417384646, 0.9221639367190003880974, 0.9466416909956290617847,
  0.9666083103968946047364, 0.9819687150345405682393, 0.9926499984472037417486,
  0.9986036451819366381565
};
static constexpr Real abscissa_n46[] = {
  0.0337721900160520415195, 0.1011624753055842395158, 0.1680911794671035286067,
  0.2342529222062697686260, 0.2993458227018700154834, 0.3630728770209957101237,
  0.4251433132828283973221, 0.4852739183881646627723, 0.5431903302618026352709,
  0.5986282897127151531771, 0.6513348462019977151064, 0.7010695120204056975121,
  0.7476053596156660540003, 0.7907300570752742551891, 0.8302468370660660530323,
  0.8659753948668580629158, 0.8977527115339419657013, 0.9254337988067539509774,
  0.9488923634460897956222, 0.9680213918539919427377, 0.9827336698041668634779,
  0.9929623489061743640730, 0.9986630421338179811282
};
static constexpr Real abscissa_n47[] = {
  0.0, 0.0660869239163556751604, 0.1318848665545148970538, 0.1971061102791118079614,
  0.2614654592149745703069, 0.3246814863377359022107, 0.3864777640846671395833,
  0.4465840731048557027251, 0.5047375838635779197742, 0.5606840059346641944830,
  0.6141786999563736085954, 0.6649877473903327291367, 0.7128889734090643016618,
  0.7576729184454386335739, 0.7991437541677419429155, 0.8371201398999021212775,
  0.8714360157968963169409, 0.9019413294385253568671, 0.9285026930123606481971,
  0.9510039692577084425898, 0.9693467873265644971457, 0.9834510030716237087649,
  0.9932552109877686346921, 0.9987187285842121091839
};
static constexpr Real abscissa_n48[] = {
  0.0323801709628693620333, 0.0970046992094626989300, 0.1612223560688917180564,
  0.2247637903946890612248, 0.2873624873554555767358, 0.3487558862921607381598,
  0.4086864819907167299162, 0.4669029047509584045449, 0.5231609747222330336782,
  0.5772247260839727038178, 0.6288673967765136239951, 0.6778723796326639052118,
  0.7240341309238146546744, 0.7671590325157403392538, 0.8070662040294426270825,
  0.8435882616243935307110, 0.8765720202742478859056, 0.9058791367155696728220,
  0.9313866907065543331141, 0.9529877031604308607229, 0.9705915925462472504614,
  0.9841245837228268577445, 0.9935301722663507575479, 0.9987710072524261186005
};
static constexpr Real abscissa_n49[] = {
  0.0, 0.0634206849826867860288, 0.1265859972696720510679, 0.1892415924618135864853,
  0.2511351786125772735071, 0.3120175321197487622078, 0.3716435012622848888637,
  0.4297729933415765246585, 0.4861719414524920421769, 0.5406132469917260665582,
  0.5928776941089007124558, 0.6427548324192376640568, 0.6900438244251321135047,
  0.7345542542374026962136, 0.7761068943454466350181, 0.8145344273598554315395,
  0.8496821198441657010348, 0.8814084455730089100370, 0.9095856558280732852130,
  0.9341002947558101490589, 0.9548536586741372335552, 0.9717622009015553801399,
  0.9847578959142130043592, 0.9937886619441677907601, 0.9988201506066353793618
};
static constexpr Real abscissa_n50[] = {
  0.0310983383271888761123, 0.0931747015600861408544, 0.1548905899981459020716,
  0.2160072368760417568472, 0.2762881937795319903276, 0.3355002454194373568369,
  0.3934143118975651273942, 0.4498063349740387891471, 0.5044581449074642016514,
  0.5571583045146500543155, 0.6077029271849502391803, 0.6558964656854393607816,
  0.7015524687068222510895, 0.7444943022260685382605, 0.7845558329003992639053,
  0.8215820708593359483562, 0.8554297694299460846113, 0.8859679795236130486375,
  0.9130785566557918930897, 0.9366566189448779337808, 0.9566109552428079429977,
  0.9728643851066920737133, 0.9853540840480058823090, 0.9940319694320907125851,
  0.9988664044200710501854
};
static constexpr Real abscissa_n51[] = {
  0.0, 0.0609611001505787247341, 0.1216954210188887669638, 0.1819770269570775453239,
  0.2415816664477987038467, 0.3002876063353319395302, 0.3578764566884095097752,
  0.4141339832263038779368, 0.4688509042860410636104, 0.5218236693661858425140,
  0.5728552163513038365223, 0.6217557046007232737550, 0.6683432211753700868644,
  0.7124444575770366445805, 0.7538953544853755257639, 0.7925417120993812052344,
  0.8282397638230648328548, 0.8608567111822923714734, 0.8902712180295273032777,
  0.9163738623097802308235, 0.9390675440029623834353, 0.9582678486139081945577,
  0.9739033680193238672317, 0.9859159917359029965838, 0.9942612604367525746210,
  0.9989099908489034951689
};
static constexpr Real abscissa_n52[] = {
  0.0299141097973387660436, 0.0896352446489005654888, 0.1490355086069491804886,
  0.2079022641563660596864, 0.2660247836050018274729, 0.3231950034348078255010,
  0.3792082691160936692466, 0.4338640677187616703090, 0.4869667456980960777824,
  0.5383262092858274383759, 0.5877586049795790699020, 0.6350869776952459242981,
  0.6801419042271677020922, 0.7227620997499831936768, 0.7627949951937449602792,
  0.8000972834304683243347, 0.8345354323267345349617, 0.8659861628460675852443,
  0.8943368905344953225208, 0.9194861289164245398935, 0.9413438536413590568435,
  0.9598318269330865525320, 0.9748838842217445031407, 0.9864461956515498406453,
  0.9944775909292160292450, 0.9989511111039502780909
};
static constexpr Real abscissa_n53[] = {
  0.0, 0.0586850543002594650227, 0.1171678090719551501402, 0.1752466621553257507297,
  0.2327214037242725936433, 0.2893939064516262064265, 0.3450688084957223566943,
  0.3995541869539529773925, 0.4526622194618457913831, 0.5042098316571334370391,
  0.5540193282770678810147, 0.6019190057137693274642, 0.6477437439165100687511,
  0.6913355756013667235407, 0.7325442308075102537821, 0.7712276549255323078660,
  0.8072524984168954782204, 0.8404945765458013754301, 0.8708392975582413516017,
  0.8981820578754266259255, 0.9224286030428121282678, 0.9434953534644418790210,
  0.9613096946231363323698, 0.9758102337149845816329, 0.9869470350233715217203,
  0.9946819193080070786360, 0.9989899477763282271210
};
static constexpr Real abscissa_n54[] = {
  0.0288167481993417776562, 0.0863545182632482152854, 0.1436054273162561539470,
  0.2003792936062135697786, 0.2564875200699973000774, 0.3117437208344682288825,
  0.3659643403721911819843, 0.4189692632552045280361, 0.4705824124813822836832,
  0.5206323343859330733270, 0.5689527681952094297316, 0.6153831983311273707299,
  0.6597693876319831246922, 0.7019638897191729193855, 0.7418265388091843162857,
  0.7792249153462540215359, 0.8140347859135678354696, 0.8461405159707729494258,
  0.8754354540655689394179, 0.9018222862847015807574, 0.9252133598666514862562,
  0.9455309751649958537638, 0.9627076457859235832569, 0.9766863288579032372000,
  0.9874206373973435585521, 0.9948751170183388849188, 0.9990266668673409838510
};
static constexpr Real abscissa_n55[] = {
  0.0, 0.0565727538183367763273, 0.1129642880593292665884, 0.1689939636468732082828,
  0.2244823006478454834001, 0.2792515532008065385498, 0.3331262788900238851894,
  0.3859339007409794297558, 0.4375052600371745918084, 0.4876751581874740972082,
  0.5362828859083432967205, 0.5831727380260321029743, 0.6281945122499281400906,
  0.6712039903198263957963, 0.7120633999866378389093, 0.7506418563480219086749,
  0.7868157811276223658976, 0.8204692985593209124542, 0.8514946066171544714602,
  0.8797923224198955060675, 0.9052718007440000257819, 0.9278514247207916968162,
  0.9474588680412107418599, 0.9640313285931351987785, 0.9775157355039892088592,
  0.9878689411988891985223, 0.9950579778474118750431, 0.9990614195648185414794
};
static constexpr Real abscissa_n56[] = {
  0.0277970352872754370940, 0.0833051868224353744402, 0.1385558468103762420128,
  0.1933782386352752582401, 0.2476029094343372039729, 0.3010622538672206690530,
  0.3535910321749545209697, 0.4050268809270912781188, 0.4552108148784595789488,
  0.5039877183843817141952, 0.5512068248555346187543, 0.5967221827706633201041,
  0.6403931068070068942679, 0.6820846126944704555015, 0.7216678344501880835225,
  0.7590204227051289022024, 0.7940269228938664980300, 0.8265791321428816516721,
  0.8565764337627486354029, 0.8839261083278275407890, 0.9085436204206554908461,
  0.9303528802474963005472, 0.9492864795619626356467, 0.9652859019054901836261,
  0.9783017091402563833769, 0.9882937155401615110899, 0.9952312260810697472163,
  0.9990943438014655843531
};
static constexpr Real abscissa_n57[] = {
  0.0, 0.0546071510016468242198, 0.1090513328087878009785, 0.1631700625912642510432,
  0.2168018287961240364136, 0.2697865731618387657633, 0.3219661683953786405902,
  0.3731848900865944585521, 0.4232898814515639509600, 0.4721316095179757095880,
  0.5195643113911876063151, 0.5654464292692367590189, 0.6096410329087153654238,
  0.6520162282809768912493, 0.6924455511995177390404, 0.7308083447445233228273,
  0.7669901193594501954917, 0.8008828945472182420760, 0.8323855211504391208290,
  0.8614039832620469447222, 0.8878516788822213295130, 0.9116496785213912127284,
  0.9327269610671016961007, 0.9510206264478767419119, 0.9664760851718866791148,
  0.9790472267094687137988, 0.9886965776502220488497, 0.9953955236784303111351,
  0.9991255656252628505720
};
static constexpr Real abscissa_n58[] = {
  0.0268470123659423558033, 0.0804636302141427293098, 0.1338482505954668570223,
  0.1868469518357613213743, 0.2393069249661534544289, 0.2910769143111091895330,
  0.3420076535979952612483, 0.3919522963307531503712, 0.4407668391868395651937,
  0.4883105372167184636155, 0.5344463096488475863998, 0.5790411351302250304899,
  0.6219664352630791110339, 0.6630984453321252664330, 0.7023185711539081134797,
  0.7395137310200422678466, 0.7745766817496527452661, 0.8074063279130881410489,
  0.8379080133393733163520, 0.8659937940748074792750, 0.8915826920220301763998,
  0.9146009285643525406865, 0.9349821375882593484808, 0.9526675575188690914428,
  0.9676062025029240901533, 0.9797550146943503091078, 0.9890790082484426364999,
  0.9955514765972909026027, 0.9991552004073866064427
};
static constexpr Real abscissa_n59[] = {
  0.0, 0.0527734840883100039516, 0.1053998790163441438369, 0.1577325055878579681152,
  0.2096255033920365449233, 0.2609342373428117116106, 0.3115157008030137003183,
  0.3612289141697948099923, 0.4099353178104189667228, 0.4574991582532666902262,
  0.5037878665577179787680, 0.5486724278083963843722, 0.5920277407040301444641,
  0.6337329662388500975125, 0.6736718645049372270240, 0.7117331186771977315959,
  0.7478106452786402318892, 0.7818038898623609056326, 0.8136181072882115714350,
  0.8431646258168722014703, 0.8703610942928822609634, 0.8951317117434720853642,
  0.9174074387881552813525, 0.9371261903534538594051, 0.9542330093769510558632,
  0.9686802216817815313536, 0.9804275739567156884496, 0.9894423651337309317821,
  0.9956996403832459646872, 0.9991833539092946837564
};
static constexpr Real abscissa_n60[] = {
  0.0259597723012477985891, 0.0778093339495365694192, 0.1294491353969450031464,
  0.1807399648734254172408, 0.2315435513760293380103, 0.2817229374232616916906,
  0.3311428482684481942523, 0.3796700565767979771549, 0.4271737415830783893074,
  0.4735258417617071111081, 0.5186014000585697474178, 0.5622789007539445391782,
  0.6044405970485103634442, 0.6449728284894770678134, 0.6837663273813554372229,
  0.7207165133557303994360, 0.7557237753065856868688, 0.7886937399322640545699,
  0.8195375261621457593685, 0.8481719847859296324905, 0.8745199226468983151293,
  0.8985103108100459419377, 0.9200784761776275528566, 0.9391662761164232494954,
  0.9557222558399961073972, 0.9697017887650527337215, 0.9810672017525981856185,
  0.9897878952222217173672, 0.9958405251188381738767, 0.9992101232274360220342
};
static constexpr Real abscissa_n61[] = {
  0.0, 0.0510589067079743493668, 0.1019846065622740689572, 0.1526442402308153005295,
  0.2029056425180584992269, 0.2526376871690534958336, 0.3017106289630307126044,
  0.3499964422040668345334, 0.3973691547257566091782, 0.4437051765385316019955,
  0.4888836222622521188206, 0.5327866265029252656384, 0.5752996513508306186003,
  0.6163117851979217247096, 0.6557160320950708716991, 0.6934095908944911554991,
  0.7292941234494651096889, 0.7632760111723121971459, 0.7952665992823596491520,
  0.8251824281086599506642, 0.8529454508476634455648, 0.8784832372148810324789,
  0.9017291624740011706420, 0.9226225813829552612575, 0.9411089866813611474775,
  0.9571401519129840913720, 0.9706742588331829082474, 0.9816760112840370796851,
  0.9901167452325170509655, 0.9959745998151202342680, 0.9992355976313634717318
};
static constexpr Real abscissa_n62[] = {
  0.0251292914218206147270, 0.0753243954962343327638, 0.1253292236158968086196,
  0.1750174592490156285556, 0.2242635856041655316689, 0.2729432026967263431886,
  0.3209333415941940040741, 0.3681127750465645296630, 0.4143623237171260481295,
  0.4595651572401133952071, 0.5036070893447559559170, 0.5463768663002510958161,
  0.5877664479530873380021, 0.6276712806468851807269, 0.6659905613354794469975,
  0.7026274922222970551209, 0.7374895252831567498607, 0.7704885960554193189909,
  0.8015413461039763715385, 0.8305693336040048513459, 0.8574992315120709228184,
  0.8822630128318973630726, 0.9047981225210934657584, 0.9250476356362037552273,
  0.9429604013923285038263, 0.9584911729739270920295, 0.9716007233716518064471,
  0.9822559490972366494911, 0.9904299711892903524280, 0.9961022963162671328851,
  0.9992598593087770296984
};
static constexpr Real abscissa_n63[] = {
  0.0, 0.0494521871161596272342, 0.0987833564469452795297, 0.1478727863578719685698,
  0.1966003467915066845576, 0.2448467932459533627484, 0.2924940585862514400361,
  0.3394255419745844024688, 0.3855263942122478924776, 0.4306837987951116006620,
  0.4747872479948043999222, 0.5177288132900332481244, 0.5594034094862850132676,
  0.5997090518776252357390, 0.6385471058213653850003, 0.6758225281149860901311,
  0.7114440995848458078514, 0.7453246483178474178293, 0.7773812629903723355633,
  0.8075354957734567600514, 0.8357135543195028434718, 0.8618464823641237195396,
  0.8858703285078534262902, 0.9077263027785315580369, 0.9273609206218432054470,
  0.9447261340410098029663, 0.9597794497589419270703, 0.9724840346975700228019,
  0.9828088105937272348625, 0.9907285468921894668108, 0.9962240127779701086021,
  0.9992829840291237803789
};
static constexpr Real abscissa_n64[] = {
  0.0243502926634244325089, 0.0729931217877990394495, 0.1214628192961205544703,
  0.1696444204239928180373, 0.2174236437400070841496, 0.2646871622087674163739,
  0.3113228719902109561575, 0.3572201583376681159504, 0.4022701579639916036957,
  0.4463660172534640879849, 0.4894031457070529574785, 0.5312794640198945456580,
  0.5718956462026340342838, 0.6111553551723932502488, 0.6489654712546573398577,
  0.6852363130542332425635, 0.7198818501716108268489, 0.7528199072605318966118,
  0.7839723589433414076102, 0.8132653151227975597419, 0.8406292962525803627516,
  0.8659993981540928197607, 0.8893154459951141058534, 0.9105221370785028057563,
  0.9295691721319395758214, 0.9464113748584028160624, 0.9610087996520537189186,
  0.9733268277899109637418, 0.9833362538846259569312, 0.9910133714767443207393,
  0.9963401167719552793469, 0.9993050417357721394569
};

struct Quadrature {
  const Real *abscissa;
  const Real *weight;
};

static const Quadrature gl_coeff[65] = {
  // n=0, n=1 are both invalid choices:
  {abscissa_n0, weight_n0}, {abscissa_n1, weight_n1},
  {abscissa_n2, weight_n2}, {abscissa_n3, weight_n3}, {abscissa_n4, weight_n4},
  {abscissa_n5, weight_n5}, {abscissa_n6, weight_n6}, {abscissa_n7, weight_n7},
  {abscissa_n8, weight_n8}, {abscissa_n9, weight_n9}, {abscissa_n10, weight_n10},
  {abscissa_n11, weight_n11}, {abscissa_n12, weight_n12}, {abscissa_n13, weight_n13},
  {abscissa_n14, weight_n14}, {abscissa_n15, weight_n15}, {abscissa_n16, weight_n16},
  {abscissa_n17, weight_n17}, {abscissa_n18, weight_n18}, {abscissa_n19, weight_n19},
  {abscissa_n20, weight_n20}, {abscissa_n21, weight_n21}, {abscissa_n22, weight_n22},
  {abscissa_n23, weight_n23}, {abscissa_n24, weight_n24}, {abscissa_n25, weight_n25},
  {abscissa_n26, weight_n26}, {abscissa_n27, weight_n27}, {abscissa_n28, weight_n28},
  {abscissa_n29, weight_n29}, {abscissa_n30, weight_n30}, {abscissa_n31, weight_n31},
  {abscissa_n32, weight_n32}, {abscissa_n33, weight_n33}, {abscissa_n34, weight_n34},
  {abscissa_n35, weight_n35}, {abscissa_n36, weight_n36}, {abscissa_n37, weight_n37},
  {abscissa_n38, weight_n38}, {abscissa_n39, weight_n39}, {abscissa_n40, weight_n40},
  {abscissa_n41, weight_n41}, {abscissa_n42, weight_n42}, {abscissa_n43, weight_n43},
  {abscissa_n44, weight_n44}, {abscissa_n45, weight_n45}, {abscissa_n46, weight_n46},
  {abscissa_n47, weight_n47}, {abscissa_n48, weight_n48}, {abscissa_n49, weight_n49},
  {abscissa_n50, weight_n50}, {abscissa_n51, weight_n51}, {abscissa_n52, weight_n52},
  {abscissa_n53, weight_n53}, {abscissa_n54, weight_n54}, {abscissa_n55, weight_n55},
  {abscissa_n56, weight_n56}, {abscissa_n57, weight_n57}, {abscissa_n58, weight_n58},
  {abscissa_n59, weight_n59}, {abscissa_n60, weight_n60}, {abscissa_n61, weight_n61},
  {abscissa_n62, weight_n62}, {abscissa_n63, weight_n63}, {abscissa_n64, weight_n64},
};

// 1D f(x1)
Real integrate(const int n, Real (*f)(Real), Real x1l, Real x1u) {
  // affine transformmation from interval [-1, 1] to [x1l, xu]
  Real m1 = (x1u - x1l)/2.0;
  Real b1 = (x1u + x1l)/2.0;
  // exploiting symmetry of weights and abscissae about x1=0.0 to only store half of coeff
  const int nelements = (n + 1)/2;
  Real sum = 0.0;
  int istart = 0;

  if (n % 2) { // odd order ---> x1=0.0 is abscissa point
    sum += gl_coeff[n].weight[0]*f(b1);
    istart = 1;  // don't double count this abscissa
  }
  for (int i=istart; i<nelements; i++) {
    Real x1 = gl_coeff[n].abscissa[i];
    sum += gl_coeff[n].weight[i]*(f(m1*x1 + b1) + f(-m1*x1 + b1));
  }
  return m1*sum;
}


// 2D f(x1, x2)
Real integrate(const int n, Real (*f)(Real, Real),
               Real x1l, Real x1u, Real x2l, Real x2u) {
  Real m1 = (x1u - x1l)/2.0;
  Real b1 = (x1u + x1l)/2.0;
  Real m2 = (x2u - x2l)/2.0;
  Real b2 = (x2u + x2l)/2.0;
  const int nelements = (n + 1)/2;
  Real sum = 0.0;
  int istart = 0;

  if (n % 2) {
    Real w_zero = gl_coeff[n].weight[0];
    sum += SQR(w_zero)*f(b1, b2);
    istart = 1;
    for (int i=istart; i<nelements; i++) {
      Real x1 = gl_coeff[n].abscissa[i];
      Real w1 = gl_coeff[n].weight[i];
      sum += w1*w_zero*(f(m1*x1 + b1, b2) + f(-m1*x1 + b1, b2));
    }
    for (int j=istart; j<nelements; j++) {
      Real x2 = gl_coeff[n].abscissa[j];
      Real w2 = gl_coeff[n].weight[j];
      sum += w2*w_zero*(f(b1, m2*x2 + b2) + f(b1, -m2*x2 + b2));
    }
  }
  for (int j=istart; j<nelements; j++) {
    Real x2 = gl_coeff[n].abscissa[j];
    Real w2 = gl_coeff[n].weight[j];
    for (int i=istart; i<nelements; i++) {
      Real x1 = gl_coeff[n].abscissa[i];
      Real w1 = gl_coeff[n].weight[i];
      sum += w1*w2*(f(m1*x1 + b1, m2*x2 + b2) + f(-m1*x1 + b1, m2*x2 + b2)
                    + f(m1*x1 + b1, -m2*x2 + b2) + f(-m1*x1 + b1, -m2*x2 + b2));
    }
  }
  return m1*m2*sum;
}


// 3D f(x1, x2, x3)
Real integrate(const int n, Real (*f)(Real, Real, Real),
               Real x1l, Real x1u, Real x2l, Real x2u, Real x3l, Real x3u) {
    Real m1 = (x1u - x1l)/2.0;
  Real b1 = (x1u + x1l)/2.0;
  Real m2 = (x2u - x2l)/2.0;
  Real b2 = (x2u + x2l)/2.0;
  Real m3 = (x3u - x3l)/2.0;
  Real b3 = (x3u + x3l)/2.0;
  const int nelements = (n + 1)/2;
  Real sum = 0.0;
  int istart = 0;

  if (n % 2) {
    Real w_zero = gl_coeff[n].weight[0];
    sum += SQR(w_zero)*w_zero*f(b1, b2, b3);
    istart = 1;
    for (int k=istart; k<nelements; k++) {
      Real x3 = gl_coeff[n].abscissa[k];
      Real w3 = gl_coeff[n].weight[k];
      for (int i=istart; i<nelements; i++) {
        Real x1 = gl_coeff[n].abscissa[i];
        Real w1 = gl_coeff[n].weight[i];
        sum += w1*w3*w_zero*(f(m1*x1 + b1, b2, m3*x3 + b3)
                             + f(-m1*x1 + b1, b2, m3*x3 + b3)
                             + f(m1*x1 + b1, b2, -m3*x3 + b3)
                             + f(-m1*x1 + b1, b2, -m3*x3 + b3));
      }
    }
    for (int k=istart; k<nelements; k++) {
      Real x3 = gl_coeff[n].abscissa[k];
      Real w3 = gl_coeff[n].weight[k];
      for (int j=istart; j<nelements; j++) {
        Real x2 = gl_coeff[n].abscissa[j];
        Real w2 = gl_coeff[n].weight[j];
        sum += w2*w3*w_zero*(f(b1, m2*x2 + b2, m3*x3 + b3)
                             + f(b1, -m2*x2 + b2, m3*x3 + b3)
                             + f(b1, m2*x2 + b2, -m3*x3 + b3)
                             + f(b1, -m2*x2 + b2, -m3*x3 + b3));
      }
    }
    for (int j=istart; j<nelements; j++) {
      Real x2 = gl_coeff[n].abscissa[j];
      Real w2 = gl_coeff[n].weight[j];
      for (int i=istart; i<nelements; i++) {
        Real x1 = gl_coeff[n].abscissa[i];
        Real w1 = gl_coeff[n].weight[i];
        sum += w2*w1*w_zero*(f(m1*x1 + b1, m2*x2 + b2, b3)
                             + f(-m1*x1 + b1, m2*x2 + b2, b3)
                             + f(m1*x1 + b1, -m2*x2 + b2, b3)
                             + f(-m1*x1 + b1, -m2*x2 + b2, b3));
      }
    }
  }
  for (int k=istart; k<nelements; k++) {
    Real x3 = gl_coeff[n].abscissa[k];
    Real w3 = gl_coeff[n].weight[k];
    for (int j=istart; j<nelements; j++) {
      Real x2 = gl_coeff[n].abscissa[j];
      Real w2 = gl_coeff[n].weight[j];
      for (int i=istart; i<nelements; i++) {
        Real x1 = gl_coeff[n].abscissa[i];
        Real w1 = gl_coeff[n].weight[i];
        sum += w1*w2*w3*(
            f(m1*x1 + b1, m2*x2 + b2, m3*x3 + b3)
            + f(-m1*x1 + b1, m2*x2 + b2, m3*x3 + b3)
            + f(m1*x1 + b1, -m2*x2 + b2, m3*x3 + b3)
            + f(-m1*x1 + b1, -m2*x2 + b2, m3*x3 + b3)
            + f(m1*x1 + b1, m2*x2 + b2, -m3*x3 + b3)
            + f(-m1*x1 + b1, m2*x2 + b2, -m3*x3 + b3)
            + f(m1*x1 + b1, -m2*x2 + b2, -m3*x3 + b3)
            + f(-m1*x1 + b1, -m2*x2 + b2, -m3*x3 + b3));
      }
    }
  }
  return m1*m2*m3*sum;
}
} // namespace GaussLegendre
